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O ■ ABSTRACT 

o 

<N . 

We discuss the space-and-time-dependent Monte Carlo code we have developed to 
^ \ simulate the relativistic radiation output from compact astrophysical objects, coupled 

£f) ' to a Fokker-Planck code to determine the self-consistent lepton populations. We have 

applied this code to model the emission from a magnetized neutron star accretion shell 
near the Alfven radius, reprocessing the radiation from the neutron sar surface. We 
explore the parameter space defined by the accretion rate, stellar surface field and the 
, level of wave turbulence in the shell. Our results are relevant to the emission from atoll 

sources, soft-X-ray transient X-ray binaries containing weakly magnetized neutron 
stars, and to recently suggested models of accretion-powered emission from anomalous 
X-ray pulsars. 

8 ■ 1. Introduction 

£3 ■ 

The high energy radiation from compact astrophysical objects is emitted by relativistic 
or semi-relativistic thermal and nonthermal leptons (electrons and pairs) via synchrotron, 
^ ■ bremsstrahlung, and Compton processes, plus bound-bound and bound-free transitions of high-Z 

elements. Since Compton scattering is a dominant radiation mechanism in this regime, the most 
efficient and accurate method to model the transport of high energy radiation is the Monte 
Carlo (MC) technique. During the past decade we have developed a versatile state-of-the-art 
space-and-time-dependent MC code to model the radiative output of compact objects (see e.g. 
Liang et al. 2000). Recently we have added the self-consistent evolution of the leptons using a 
Fokker-Planck scheme. The lepton evolution is then coupled to the photon transport. Since this 
is the first time that we report on results obtained with this code, the first part of the present 
paper is devoted to a detailed description of the capabilities of the code (§2) and its verification in 
comparison with previous work (§3). 
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In the second part of this paper (§4), we apply this coupled MC-FP code to the study of 
the reprocessing of blackbody radiation from a magnetized neutron star surface by a magnetized 
spherical shell at the Alfven radius. Recent observations of weak-field neutron star binaries, such 
as low-luminosity, X-ray bursters (atoll sources, e.g., 4U 1608-522: [Zhang et al. 1996| , 4U 1705-44: 
Barret et al. 1996| ), bursting soft X-ray transients (e.g., Aql X-l: Harmon et al. 1996] ), or pulsar 



binary systems (e.g., PSR B1259-63: Tavani et al. 1996|) indicate that many of them exhibit soft 



(photon index ;> 2) power law tails extending beyond ~ 100 keV, at least episodically, in addition 
to the thermal component at temperatures of ~ a few keV, which presumably originates from the 
stellar surface. The luminosity of this high-energy tail appears to be anti-correlated with the soft 



X-ray luminosity (Barret & Vedrenne 1994, Tavani & Liang 1996) 



The origin of this high energy tail is unexplained at present. It could be due to thermal 
Comptonization by a hot coronal plasma, or it could be due to nonthermal emission. Tavani and 
Liang (1996) examined systematically the possible sites of nonthermal emissions and concluded 
that the Alfven surface is the most likely candidate since the dissipation of the rotation energy of 
the disk is strongest there, due to magnetic reconnection and wave turbulence generation. Here 
we first focus on particle acceleration by wave turbulence. 

Because of the 1-D nature of our code at this stage, the neutron star is assumed to emit 
isotropic blackbody radiation, the reprocessing magnetospheric plasma is idealized as part of 
a spherical shell at the Alfven radius, and the magnetic field is taken to be nondirectional, so 
that the radiation output is isotropic. Even though this is not a perfect representation of a 
quasi-dipolar magnetosphere and accretion flow, we believe that, except for very special viewing 
angles, our output spectra should be a reasonable approximation to the angle-averaged output of 
the reprocessing by a hybrid thermal/nonthermal magnetosphere. We assume that the leptons 
are energized by Coulomb collisions with virial ions and accelerated nonthermally by Alfven and 
whistler wave turbulence, and cooled by cyclotron/synchrotron, bremsstrahlung, and inverse 
Comptonization of both internal soft photons and blackbody photons from the stellar surface. 

The primary focus of the parameter study presented in §4 is the application to weakly 
magnetized neutron stars with surface magnetic fields of -B S urf Ss 10 11 G. In particular, we will 
show that the anti-correlation of the hardness and luminosity of the hard X-ray emission with 
the soft X-ray luminosity is a natural consequence of the energetics of particle acceleration and 
cooling near the Alfven radius. We predict that the nonthermal tails in the hard X-ray spectra of 
accreting, weakly magnetized neutron stars may extend up to ~ 1 MeV. A solid detection and the 
measurement of the cutoff energy of these high-energy tails by the INTEGRAL mission, scheduled 
for launch in 2002, will provide important constraints on accretion-based models for the hard 
X-ray emission from accreting neutron stars. 

In this context, it is interesting to note that Chatterjee, Hernquist, & Narayan (2000; see also 
Mcrcghctti & Stella 1995, Wang 1997, |Chattcrjcc fc Hernquist 2000 ) have recently proposed a 



similar type of accretion-powered emission for anomalous X-ray pulsars (AXPs), as an alternative 



-3 - 



to models based on magnetic- field decay ( [Thompson &z Duncan 1996 ) or residual thermal energy 
QHeyl fc Hernquist 19971) . According to Chatterjee et al. (2000) the X -ray emission from AXPs 
(which generally consists of a soft, thermal component with kT ~ 0.3 - 0.4 keV plus a hard 
X-ray tail with photon index T ~ 3 - 4) is powered by accretion of material from the debris of 
the supernova which had formed the neutron star, onto a the surface of the neutron star, which 
possesses a typical pulsar magnetic field of -B sur f ~ 10 12 G. Therefore, we extend our parameter 
study to parameter values relevant to accreting pulsars. However, we point out that in the case of 
a magnetic field as high as -B SU rf ~ 10 12 G, the assumed shell geometry and the quasi-isotropy of 
the emission from the neutron star surface may be a gross over-simplification. However, although 
consequently the precise parameter values used in this region of the parameter space should 
not be taken at face value, our parameter study might still provide interesting insight into the 
dependence of the equilibrium electron and photon spectra on the various input parameters in the 
high-magnetic-field case. 



2. Physics of photon and lepton evolution 

We use the Monte Carlo (MC) technique (Podznyakov, Sobol, &: Sunyaev 1983; Canfield, 



Howard, & Liang 1987; Liang 1993; Hua, Kazanas, & Titarchuk 1997) to simulate relativistic 
photon transport. We include the full (energy-and-angle-dependent) Klein-Nishina cross section 
for Comptonization, relativistic bremsstrahlung from lepton-ion and lepton-lepton scattering 



( Dermer 1984 ), and cyclo-synchrotron processes ( Brainerd 1984 ) in magnetic fields. The MC 
photon transport is fully space-and-time-dependent. Photons are born with a certain "weight" 
(which is basically the total energy in photons within a spatial region or emitted from a photon 
source, divided by the number of Monte-Carlo particles) which is diminished by absorption and 
escape, until the weight drops below a user-specified limit, at which point the photon is "killed". 
In the simulations shown later in the paper, we specify this energy weight cutoff as 1/100 of 
the initial statistical weight of the photon. The final results are insensitive to the choice of this 
energy weight cutoff, as long as it is <C 1. Surviving photons are sampled at boundaries to provide 
time-and-frequency-dependent spectral output. In addition to self-emitted photons from the 
plasma, soft photons can be injected at zone boundaries and inside volume elements. Currently 
the code can handle 1-D spherical, cylindrical or slab geometries with an arbitrary number of 
spatial zones. However since the photon ray tracing is done with full angle informations, the 
generalization to 2-and-3-D transport is straight forward. The maximum number of photon 
frequency bins is 128. For more details of this code see Canfield et al. (1987) and Bottcher & 
Liang (1998). A typical MC run with a million particles at a Thomson depth of a few takes 
10s of minutes on a DEC alpha server. Since the CPU time usage scales as the square of the 
Thomson depth ( number of scatterings), Thomson thick runs can be quite time consuming. We 
are currently developing a random walk scheme for Wien photons trapped in Thomson thick 
zones, which would save large amounts of CPU time without introducing too much error. 
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The lepton population is computed locally in each spatial zone using the Fokker-Planck 
approximation (Dermer, Miller, & Li 1996; Li, Kusunose & Liang 1996a), taking into account 
coulomb and Moller scattering, stochastic acceleration by Alfven and whistler wave turbulence, 
and radiative cooling (plus pair processes if necessary). In general the lepton population consists 
of a low energy thermal population plus a high-energy tail truncated at the highest energies by 
radiative cooling. The photon and lepton evolutions are coupled to each other via a quasi-implicit 
time scheme in which we use an average of the photon-electron energy exchange rates between 
two subsequent time steps to compute the electron cooling rates, in particular due to Compton 
scattering. The Fokker-Planck equation governing the electron evolution is solved using a fully 
implicit scheme. In contrast to several other codes currently used in the literature (e.g., [Stern et 



aL 1995| ; |Li et al. 1996al pVlalzac fc Jourdain 2000| ), we solve the evolution of the entire electron 



population with the FP scheme and do not introduce any artificial separation between thermal and 
non-thermal particles. Since the lepton distribution typically evolves much faster than the photon 
distribution, each photon cycle contains many lepton cycles. The user, however, can always turn 
off the nonthermal lepton acceleration and assume a strict thermal population whose temperature 
can be computed self-consistently from energy balance alone. 

In oder to calculate the emissivities and opacities for thermal cyclotron, non-thermal 
synchrotron, and thermal bremsstrahlung emission (and only for this purpose), the electron 
distribution calculated with our Fokker-Planck scheme is decomposed into a thermal population 
plus a non-thermal tail. For the thermal bremsstrahlung and non-thermal synchrotron emission 
and absorption, the standard expressions from Rybicki &; Lightman (1979) are used. For the 
thermal cyclotron emission, we add explicitly over the first 5 harmonics, beyond which we use the 
asymptotic continuum representation of Mahadevan, Narayan & Yi (1996). 

Since in many cases of interest to the current investigation the moderately to strongly 
magnetized coronal plasma is optically thick at low frequencies, hv <C 1 keV, due to synchrotron- 
(self-)absorption, we use the following simplification in order to save CPU time: For any given 
frequency hu <^ 1 keV, we compare the absorption length to the Compton mean free path, 
^Compt an d t ne radial extent AR of the current zone of the Comptonizing region. If 

^ bs <0.1min{Z£ ompt) Ar] (1) 

any soft photon produced at this frequency has a very small probability of escaping the current 
zone or being up-scattered (to a frequency at which the absorption length will be different so that 
the above criterion has to be re-evaluated after the scattering event) before being re-absorbed. 
Thus, in the frequency ranges where Eq. ([!]) is fullfilled, the radiation escaping at the boundary 
of that zone will be approximately given by the respective section of the thermal blackbody 
spectrum. Accordingly, our code neglects the volume emissivities in those frequency ranges and, 
instead, produces thermal blackbody photons at the zone boundaries. 



During each photon time step, the code keeps track of the energy transferred between photons 
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and electrons due to Compton scattering and cyclotron/synchrotron and bremsstrahlung emission 
and absorption. The respective heating and cooling rates are used to scale the energy loss/gain 
coefficients of electrons at a given energy E e = "fm e c 2 . For the bremsstrahlung energy loss rate, we 
use the approximate scaling law (dj/dt)^ oc — 7 L1 (Bottcher, Pohl, & Schlickeiser 1999). Coulomb 
heating/cooling is included using the energy exchange and dispersion rates given in Dermer & 
Liang (1989). We assume that the ions have a pre-specified temperature kT p ^ 1 MeV, which is 
not significantly affected by any changes of the electron temperature. For the energy exchange 
due to M0ller scattering (elastic electron-electron scattering) we use the electron energy exchange 
and dispersion coefficients given in Nayakshin & Melia (1998). 

In addition to Coulomb interactions and radiative losses, we also account for stochastic (2 nd 
order Fermi) acceleration due to Alfven and Whistler waves. For a given background magnetic 
field B, the magnitude and spectrum of hydromagnetic plasma wave turbulences are determined by 
the parameters 5 2 = (AB/B) 2 , where AB is the amplitude of the magnetic-field fluctuations, and 
q, the spectral index of the turbulence spectrum. We will generally use q = 5/3, characteristic of 
Kolmogorov turbulence. The electron acceleration and energy dispersion rates are calculated using 
the formalism of Schlickeiser (1989). However, we have to take into account that those plasma 
waves interacting efficiently with the low-energy, quasi-thermal part of the electron spectrum will 
be efficiently damped and absorbed by the process of Landau damping (e.g., Schlickeiser, Fichtner 
& Kneller 1996). This leads to a strong truncation of the Kolmogorov wave spectrum above a 
critical wave number for which the equivalent absorption depth through the region occupied by 
the hot plasma exceeds unity. We take the effect of Landau damping into account by introducing 
an effective absorption depth at wave number k, and correcting the acceleration rate 7^ of 
Schlickeiser (1989) due to the "optically thin" plasma wave spectrum by an absorption term: 



n (1 — e Tfc res) 

1A=7°A K — (2) 

&res 

Here r / t rcs = Th rcs t\, and is the Landau damping rate at wave number fc rcs at which Alfven 
waves are resonating preferentially with electrons of energy 7. = AR/va is the Alfven crossing 
time of the zone. This modified acceleration term can be renormalized either to correspond to 
the pre-specified value of S 2 for electrons resonating with waves in the weakly damped part of the 
Alfven wave spectrum (i.e. at high electron energies), or by specifying a heat input rate into the 
electron ensemble due to resonant wave/particle interactions. 

As an option, our code can account for pair production and annihilation and the corresponding 
photon absorption and emission. The pair annihilation rates and annihilation radiation emissivities 
are taken from Svensson (1982), and the 77 pair production rate of Bottcher & Schlickeiser 
(1997) is used. In addition, simple Compton reflection schemes, using the Green's functions for 
reflection off neutral disk material of White, Lightman & Zdziarski (1988) and Lightman & White 
(1988), can be used to simulate a Compton reflection component either off the inner boundary 
(corresponding to a quasi-homogeneous slab geometry) or reflecting part of the radiation escaping 
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at the outer boundary (e.g., corresponding to a cold outer disk). 

Assuming local isotropy of the electron distribution, we solve the one-dimensional 
Fokker-Planck equation 

+ \^[n e {l,t)D{ 1 ,t)\, (3) 

where D is the energy dispersion coefficient. To solve Eq. |3], we use an implicit version of the 
discretization scheme proposed by Nayakshin & Melia (1998). We choose a logarithmic spacing in 
electron kinetic energy, x, = 7, — 1. In the following the subscript i refers to electron energy, while 
the superscript n refers to time. We define /" = n e (^i,t n )/n e The discritization is then given by 



dt 



d 



07 
~dt 



aAx^Ax^i + Axi) ) 

where a; = (dj/dt)(xi), Di = [d(A'j) 2 /dt](xi), a = 2/(1 + q), p = 2q/(l + q), with q = x i+ i/xi. 
The system of equations is supplemented by the boundary conditions specified in Appendix A of 
Nayakshin & Melia (1998). The system of equations (|j) is in tridiagonal form and can be easily 
solved to find the electron distribution at time t n+l . In order to evaluate the energy exchange and 
dispersion coefficients a™ +1 and -D™ +1 the code performs an energy balance calculation to find the 
average electron energy in the ensemble after the current time step, which is then used to calculate 
the energy exchange and dispersion coefficients due to M0ller scattering. All other coefficients 
evolve sufficiently slowly (i.e. on time scales much longer than the electron-evolution time scale) 
so that the coefficients evaluated under the conditions at the beginning of the current time step 
can be used. 

The implicit scheme of Eq. |I] is known to approach the equilibrium solution for the electron 
distribution exactly, although the temporal evolution calculated this way becomes inaccurate 
on short time scales. Since in the problems of interest here, the electron distributions typically 
evolve on timescales much shorter than the photon distributions, the electrons always attain a 
distribution close to local equilibrium during each photon time step. Hence, the degree of accuracy 
provided by the implicit scheme is sufficient for our purposes. As an option, the code can perform 
the same simulations using an explicit scheme to solve the Fokker-Planck equation. In all test 
cases, both schemes produced virtually identical results, but the implicit scheme executes a factor 
of ~ 10 - 100 faster because of the larger individual Fokker-Planck time steps allowed in this 
scheme. 

However, we did not find an appropriate method to calculate the coefficients for the pair 
production and annihilation rates at time step n + 1. For this reason, we have to correct the 
electron and positron distribution functions for pair production and annihilation in an explicit 
manner. This becomes obviously inaccurate in the case of strongly pair producing model 



/ „n+l f n+l _ n+l rn+1 \ / 

I" = /f +1 + At ( +1 Ai'+'ai".', )~ At { 
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situations. Thus, our scheme is a significant improvement over existing schemes only for pair 
deficient situations with pair fractions / pa i r = n e +/n p <; a few %. For strongly pair dominated 
situations, we would have to adjust the individual time steps of our simulations so far that such 
simulations would take a prohibitive amount of computing time. 

In the present paper, we focus on results for equilibrium situations, for which we let the 
code evolve until a stable electron distribution and photon output spectrum is reached. The 
convergence of the electron and photon spectra can be greatly accelerated if one starts out with 
an appropriate first guess for the equilibrium electron temperature. For this reason, we have 
developed an analytical estimate of the equilibrium electron temperature, which is described in 
Appendix A. This analytical estimate is implemented in our code and can be used to determine 
the appropriate initial conditions for our equilibrium simulations. Detailed tests and applications 
of the time-dependent features of the code will be presented in future publications. 



3. Tests of the numerical scheme and comparison with Previous Results 

For verification of our code, we have first compared the individual energy exchange and 
diffusion rates, drf/dt and D(j), with those obtained in earlier work. Our numerical energy 
exchange and diffusion rates due to Coulomb scattering are in agreement with those of Dermer 
& Liang (1998)Q The numerical values of the M0ller scattering energy exchange and diffusion 
coefficients were found in good agreement with those of Nayakshin & Melia (1998)Q. 

We have verified that in cases in which radiative cooling and hydromagnetic acceleration 
are inefficient, our Fokker-Planck scheme correctly produces a thermal electron distribution in 
temperature equilibrium with the protons. 

In the analysis of coronal energy and radiation transfer, it is customary to paramatrize 
energy input and dissipation rates u [ergs cm -3 s" 1 ] by the dimensionless compactness, I, which in 
spherical geometry, is defined by l sp h = Atktt R 2 «/(3m e c 3 ), where R is the radius of the spherical 
volume. In slab geometry, this becomes Z s i a b = <tt H 2 ii/(m e c 3 ), where H is the thickness of the 
slab. In our simulations, we specify the heating mechanisms to be Coulomb heating through the 
thermal protons and resonant wave-particle interactions. Thus, leaving the shape of the electron 
spectrum general, we calculate the respective energy dissipation rates as 



""Coiil/ A — 




3 note that in their Eq. (A2), the last 7* has to be replaced by (7*) 2 

4 note that there is a '-' sign missing in front of the term oc (7 — 71) 2 in their Eq. (35) 
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where the subscript 'Coul/wp' stands for Coulomb scattering and wave-particle interaction, 
respectively. Obviously, the resulting compactness will depend not only on the energy density in 
protons and magnetic turbulence, but also on the current electron distribution (i.e. temperature, 
if the distribution is predominantly thermal), and can thus not be specified a priori as a free 
parameter. For comparison with previous work, we need to find appropriate values of the 
ion temperature and the wave turbulence amplitude in order to reproduce the dissipation 
compactness quoted in those papers. We define the temperature corresponding to the electron 
spectrum resulting from our simulations by requiring that the average electron Lorentz factor, 
(7) = (l/n e ) / d7 7n e (7), be equal to the average Lorentz factor of a thermal electron population 
of the respective temperature, 



To test our energy balance calculations, we have performed a series of simulations to reproduce 
the results of the slab corona model of Dove et al. (1997), in particular the inset of their Fig. 
2, with which we find good agreement. For those simulations, we specified ion temperatures of 
80 MeV <^ kTi <J 250 MeV, a magnetic field in equipartition with the ion population, and a 
negligibly low level of Alfven wave turbulence to reproduce various values of the local heating 
compactness l c . As mentioned earlier, our code is rather inefficient in simulating strongly 
pair-producing, high-Z c situations because of the explicit scheme to solve for pair balance. For this 
reason, we restrict the applications of the code in its present version to model situations with pair 
fractions <^ a few %. 

Li et al. (1996a, b) have developed a code solving simultaneously the Fokker-Planck equations 
for both the electron and the photon distribtions in a homogeneous medium, including Coulomb 
interactions, thermal bremsstrahlung, Compton scattering, resonant wave-particle interactions, 
and pair production and annihilation. There are two major differences between their approach 
and ours: (1) We use a Monte-Carlo method to solve the photon transport, and (2) we solve 
the Fokker-Planck equation for the entire electron spectrum, while Li et al. (1996a, b) split the 
electron distribution up into a thermal "bath" plus a non-thermal tail, assuming a priori that 
electrons of energies 7 < 7 t h r = 1 + 40 e attain a thermal distribution, and that acceleration due to 
wave/particle interactions affects only particles beyond 7thr- The latter simplification is justified 
by the argument that at low electron energies the thermalization time scale is much shorter than 
any other relevant time scale, and that long-wavelength plasma waves, with which low-energy 
electrons resonate preferentially, are strongly damped and in energy elequilibrium with the thermal 
pool of electrons. In our simulations, both effects are taken into account self-consistently without 
making the a-priori assumption of the existence or development of a non-thermal population. 

Li et al. (1996b) present two model calculations to explain hard tails observed in the 
hard X-ray / soft 7-ray spectra of Cyg X-l and GRO J0422+32. A spherical region of radius 




(6) 
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R = 1.2 x 10 8 cm is assumed. For the case of Cyg X-l, they specify a total heating compactness of 
I = 4.5, and a non-thermal heat input into suprathermal particles (7 > 7thr) of l s t/l = 0.15, wich a 
turbulence amplitude of 5 2 = 0.059. A soft blackbody radiation component at kT s = 0.1 keV from 
the outer boundary of the sphere is assumed to provide a soft radiation compactness of l s /l = 0.07. 
The seed Thomson depth of the region is r p = 0.7. In their simulations, Li et al. (1996b) find a 
significant suprathermal tail in the resulting electron distribution, which leads to an excess hard 
X-ray / soft 7-ray emission, consistent with the observed one. The temperature of the thermal 
part of their electron ensemble is found at kT e = 139 keV, and the equilibrium pair fraction is 
/pair ps 1.8 %. For the case of GRO J0422+32 they specify I = 0, l st /l = 0.06, and l s /l = 0.08. 
This resulted in a weaker nonthermal electron tail, an equilibrium temperature of kT e = 133 keV, 
a wave amplitude of 5 2 = 0.065, and a pair fraction of / pa j r ss 3 %. 

For comparison with our code, we have run simulations with the same total compactness I, 
soft compactness l a , soft blackbody temperature kT^B, radius R, seed Thomson depth r p , and 
the same values of the plasma wave amplitude normalization 5 2 . A difficulty in comparing the 
two codes is that in Li et al. (1996a, b) the ion temperature and magnetic field are not specified 
explicitly. In our comparative simulations we assume a magnetic field in equipartition with the 
ion energy density. The results of our simulations are illustrated in Figs. [l| and |2[ Our results 
are qualitatively similar to those of Li et al. (1996a, b). However, we find somewhat lower 
equilibrium temperatures and stronger nonthermal electron tails as well as stronger suprathermal 
acceleration compactnesses l s t- The lower temperatures may be attributed to a rather prominent 
cyclotron/synchrotron cooling (comparable to Compton cooling) in our simulations (see Figs, [ljb 
and|2]b). This, in combination with the larger Z st values, seems to indicate that in our simulations 
we have used higher magnetic field values than Li et al. (1996a, b). However, even abandoning 
the equipartition assumption, we could not find self-consistent parameter values resulting in the 
same combination of input parameters used by Li et al. (1996a, b). Given this descrepancy in 
the way of input parameter specification, the qualitative agreement between the results of the two 
codes, using very different numerical methods, is encouraging. 



4. Models of accretion onto a magnetized neutron star 

While in the first part of this paper we were describing the general features and the verification 
of our MC/FP code, we are now applying this code to model the electron dynamics and photon 
transport arising from models of accretion onto a magnetized neutron star. In both weakly 
magnetized neutron stars (atoll sources and soft X-ray transient neutron star binary systems) 
with -Bgurf <^ 10 11 G and X-ray pulsars (including, possibly, anomalous X-ray pulsars) with 
-Ssurf ~ 10 12 G, energy dissipation will be most efficient at the Alfven radius, where the optically 
thick, geometrically thin outer accretion disk is disrupted and the dynamics of the accretion flow 
becomes dominated by the magnetic field. We idealize this region of efficient energy dissipation 
at the Alfven radius of disk accretion onto a magnetized neutron star as part of a spherical shell 
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whose distance tq, magnetic field, and column thickness are fixed by the accretion rate (Ghosh & 
Lamb 1979a] , b, 1990). The distance r$ is given by 



r = 2 x 10 8 / fi% 7 l* 2/7 M~ 1/7 R 6 2/7 cm (7) 

where //3o = (neutron star magnetic moment)/(10 30 Gem 3 ), = L/L^dd, = M^q/Mq, and 

Rq = i?Ns/(10 6 cm). For the current simulations, for definiteness, we fix the Ghosh-Lamb fudge 
parameter / = 0.3, and set M* = Rq = 1. Hence, the dipole magnetic field at the Alfven radius is 

B = 4.2 x 10 6 lT Msq 577 Ml 11 Rf f Q l G, (8) 
where /0.3 = //0.3. The virial ion temperature at tq is 

m = V^^L « 9.3 f^V 1 M* MeV. (9) 



3 ro \10 7 cm 
The column density of the shell can be estimated using the poloidal accretion rate 
M ~ Atttq Aro Hi v p run, where we assume that the poloidal velocity v p ~ vg/2 with vg 
being the free-fall velocity. Hence, the radial Thomson depth of the shell is approximately: 

A Ma T ^ nnv r8/7 -2/7 , ,4/7 D l/7 ,-1/2 / ln x 

r T = Ar rij cr T ~ psO.97// /i 30 MJ R 6 ' f 03 < . (10) 

The neutron star (taken to be a 10 km spherical surface) is assumed to emit a blackbody luminosity 
at the temperature /cTbb = 1.78 Z* keV. In addition the level of wave turbulence is specified by 
the dimensionless amplitude 5 2 = (AB/B) 2 and the spectral index q. The minimum wavevector 
k m i n is set to 2-7r/(Aro), where Aro ~ O.Itq is the shell thickness ( phosh &, Lamb 1979a , b, 1990). 



Due to spherical symmetry of our simulations the magntic field is assumed to be nondirectional in 
the shell and the synchrotron emissivities and absorption coefficients are angle- averaged. 

We have explored the parameter space by varying the accretion rate — corresponding to a 
variation of the parameter I* — , the magnetic field — corresponding to a variation of /U30 — , and 
the amplitude of wave turbulence, 5 2 . q is set = 5/3 in all runs presented in this paper. 

In Fig. |, we demonstrate the effect of a varying accretion rate in the case of no turbulence, 
5 2 = 0, and for a fixed magnetic moment of the neutron star, fi^Q = 1, corresponding to a strong 
surface field of -B S urf ~ 10 12 G. As the accretion rate decreases, the Alfven radius moves outward, 
implying a lower magnetic field at the Alfven radius, and a lower ion temperature and Thomson 
depth of the active shell. At the same time, however, Compton cooling becomes less efficient 
due to the reduced soft photon luminosity of the neutron star surface and due to the larger 
distance of the active region from the surface. This reduction of the soft photon compactness 
leads to an increasing equilibrium electron temperature in the accretion shell. Consequently, as 
the accretion rate decreases, the photon spectra change in the following way: For /* = 1, the 
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hard X-ray spectrum smoothly connects to the peak of the thermal blackbody bump from the 
neutron star surface, and shows a quasi-exponential cutoff at high energies. For lower Z*, the 
spectrum turns into a typical low-hard state spectrum of X-ray binaries with the thermal bump 
clearly distinguished from a hard power-law + exponential cutoff at high energies. The hard X-ray 
power-law becomes harder with decreasing Z*. Fig. || shows the same sequence of decreasing I* 
for = 10~ 3 , corresponding to -B S urf ~ 10 9 G. At high accretion rates, a strong reduction of 
the equilibrium electron temperature with respect to the strong-field case results, leading to a 
softening of the photon spectrum with decreasing This is due to the increasing magnetic field 
at the Alfven radius as the neutron star magnetic moment decreases (because the Alfven radius 
decreases, overcompensating for the decreasing magnetic moment, see Eq. |8|), resulting in a lower 
electron temperature caused due to increasing cyclotron/synchrotron cooling. For low accretion 
rates (7* a few %) the electron and photon spectra for the strong-field case and the low-field 
case are virtually indistinguishable from each other. 

Figs. H and ^ illustrate the effects of a varying accretion rate and turbulence level on 
the electron and photon spectra. An increasing turbulence level, obviously, leads to a more 
pronounced nonthermal tail or bump in the electron spectrum. At the same time, as a result of 
increased Compton cooling on synchrotron photons produced by the nonthermal electrons, the 
temperature of the quasi-thermal portion of the electron spectrum decreases as the turbulence level 
is increasing. Consequently, the photon spectrum at soft to medium-energy X-rays becomes softer 
with increasing 5 2 , while at hard X-ray and 7-ray energies an increasingly hard tail develops. As 
in the case without turbulence, a decreasing accretion rate leads to a higher electron temperature 
and the transition from a smoothly connected thermal blackbody + thermal Comptonization 
power-law spectrum to a typical low/hard state X-ray binary spectrum. 

Fig. [?] demonstrates the moderate dependence, in particular of the resulting photon spectra, 
on the magnetic moment of the neutron star for an intermediate value of the accretion rate, 
I* = 0.25. Nonthermal tails in the electron spectra become more pronounced with increasing 
neutron star magnetic moment. For the weak-field case with = 1CP 3 , we find that the 
high-energy end of the electron spectra are always truncated with respect to a thermal distribution 
due to strong Compton losses, rather than developing a non-thermal tail. 

At very low accretion rates, resulting in very low proton and electron densities in the accretion 
shell, we find that even at low turbulence level the heating due to stochastic acceleration strongly 
dominates over Coulomb heating. Consequently, the resulting electron spectrum becomes strongly 
non-thermal. We need to point out that in those cases, our treatment is no longer self-consistent 
since in our code the attenuation of Alfven waves is calculated under the assumption that it is 
dominated by Landau damping in a thermal background plasma. 

The results illustrated in Figs. |3| and ^ are in excellent agreement with the general trend 
(Barret & Vedrenne 1994, Tavani & Liang 1996) that hard tails in LMXBs and soft X-ray 
transients containing weakly magnetized neutron stars are only observed during episodes of low 
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soft-X-ray luminosity. While at large soft-X-ray luminosity, the X-ray emission out to ^> 20 keV 
is dominated by the thermal component from the neutron star surface and the hard tail is very 
soft with a low cut-off energy at E c ~ 100 keV, in lower-luminosity states the hard X-ray tail 
becomes very pronounced with photon indices T ~ 2 - 3 and extends out to E c <J 1 MeV. 
Dedicated deep observations by the upcoming INTEGRAL mission should be able to detect this 
hard X-ray emission from neutron-star-binary soft X-ray transients and atoll sources in both high 
and low soft-X-ray states and thus provide a critical test of this type of accretion model for weakly 
magnetized neutron stars. 

Interestingly, a similar anti-correlation of the hard X-ray spectral hardness with the soft 
X-ray luminosity seems to exist for high surface magnetic fields (see Figs. || and ||) only in the 
case of very low accretion rates (7* ~ 0.01) and rather high magnetic turbulence levels (5 2 ;> 0.01). 
The hard X-ray spectral indices resulting from our simulations are in excellent agreement with the 
values of r ~ 3 - 4 generally observed in anomalous X-ray pulsars. This may provide additional 
support for accretion-powered emission models for AXPs. Our simulations predict cutoff energies 
of E c ~ 100 - 500 keV. 



5. Summary and Conclusions 

We have reported on the development of a new, time-dependent code for radiation 
transport and particle dynamics. The radiation transport, accounting for Compton scattering, 
bremsstrahlung emission and absorption, cyclotron and synchrotron emission and absorption, 
and pair processes, is done using a Monte-Carlo method, while the electron dynamics, including 
radiative cooling, Compton heating/cooling, and stochastic acceleration by resonant interaction 
with Alfven/whistler wave turbulence, are calculated using an implicit Fokker-Planck scheme, 
coupled to the Monte-Carlo radiation transfer code. 

In the second part of this paper, we have applied our code to the static situation of a shell at 
the Alfven radius of a magnetized neutron star. This situation is representative for accretion onto 
weakly magnetized neutron stars (bursting atoll sources or soft X-ray transients containing wekaly 
magnetized neutron stars) as well as to recently suggested models of accretion-powered emission 
from anomalous X-ray pulsars. 

The main results of our parameter study are: 

(1) The lepton thermal temperature increases, and the hard X-ray photon spectrum becomes 
harder as the accretion rate is decreasing. At the same time, the normalization of the hard X-ray 
power-law, relative to the thermal blackbody from the neutron star surface, becomes smaller. 

(2) The nonthermal tails in the electron and photon spectra become more dominant and 
harder as the turbulence level is increasing. At the same time, the quasi-thermal electron 
temperature decreases, leading to a softer hard- X-ray spectrum. 
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(3) For low accretion rates (<; a few %), the photon spectra are only very weakly dependent 
on the magnetic moment of the neutron star. For higher accretion rates, an increasing neutron 
star magnetic moment leads to a moderate hardening of the hard X-ray spectrum. 

Our results are in good agreement with the non-detection of hard X-ray tails during soft 
X-ray high states of systems containing weakly magnetized neutron stars. However, we predict 
that a hard X-ray excess beyond ~ 20 keV should exist even in the high/outburst state. This 
excess should have a cutoff energy of E c ~ 100 - 200 keV. In the low/quiescent state, the predicted 
spectral indices of T ~ 2 - 3 are in excellent agreement witht he observed hard X-ray excesses 
observed from sources believed to contain weakly magnetized neutron stars. We predict that these 
spectra should extend out to E c <^ 1 MeV. 

In the strong-field case, representative of a recently suggested model for anomalous X-ray 
pulsars, we expect a strong correlation between the hard X-ray spectral hardness and the soft 
X-ray luminosity only in the case of very low accretion rates and high magnetic turbulence level. 
The predicted hard X-ray spectral indices are generally in very good agreement with the observed 
values of T ~ 3 - 4. Cutoff energies of E c ~ 100 - 500 keV are predicted. 

Dedicated, deep observations by the upcoming INTEGRAL mission should be able to 
establish the existence of the hard power-law tails predicted in the models discussed here, and to 
constrain the high-energy cutoff of these tails. These measurements will be essential ingredients 
for a more detailed modeling of the physical conditions governing the accretion onto magnetized 
neutron stars. 

The work of MB is supported by NASA through Chandra Postdoctoral Fellowship Award No. 
9-10007, issued by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical 
Observatory for and on behalf of NASA under contract NAS 8-39073. This work was partially 
supported by NASA grant NAG5-4055. We wish to thank the referee for very helpful comments, 
and D. Marsden for pointing out the relevance of our work to anomalous X-ray pulsars. 

A. Estimate of equilibrium temperature 

The code outlined above can be used for equilibrium situations by simply letting it evolve 
until both the photon and electron distributions have relaxed to a steady state. In order to provide 
a consistency test of our numerical methods, we compute a quasi- analytical estimate for the 
expected equilibrium electron temperature. In our equilibrium simulations, we use this estimated 
equilibrium temperature as initial condition in order to accelerate the convergence. 

Assuming that the electron distribution is roughly thermal, we may estimate the heating and 
cooling rates due to the various processes as follows. Let 
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be the average energy per electrons. We approximate the cyclotron emissivity by the high-frequency 
continuum limit given in Mahadevan et al. (1996): 



where e = kT e /(m e c 2 ), v = v/{y c ® 2 ), and v c = eB/{2n m e c). 



9v\^ 
2 



(A2) 



Denoting by AR the thickness of the emitting shell or slab, the synchrotron-self-absorption 
frequency is then determined by solving for 

1 = AR KSSA = AR^\. (A3) 
The cyclotron/synchrotron cooling rate may then be estimated as 

(^SSA oo \ 

Surface / dvB u {Q e )+V J dv j„(0 e ) , (A4) 
^SSA / 

where V and ^4 S urface are the volume and the surface area of the emitting shell or slab. 

Assuming that the photon field inside the emitting volume is dominated by low-energy 
photons with e <C © e , where e = hv/(m e c 2 ), the Compton cooling rate may be approximated as 

( dw \ a ^ ^ 3 te) / * _ 

\dtJ IC k 2 (£) 

The photon energy density u p h is calculated iteratively from the contributions of the various soft 
photon fields, and repeated Compton scatterings: Assume that the soft photon energy density 
u s = u sy + ut>r + "Uext (sycnhrotron + bremsstrahlung + external radiation field) is centered around 
a soft photon energy e s , and that tt ^ 1- The average energy change of a photon with mean 
photon energy e& (Compton scattered k times) changes on average by a factor Ae k = e k (4© — £&) 
in the course of the k + 1. scattering. Then, the total internal photon energy density is 

(oo n— 1 \ 
1 + E T T II [ 40 e " ^ , (A6) 
n=l k=0 ) 

where e = e s and e fc+ i = e k (1 + 4© e - e k ). 
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In the cases we are interested in, the bremsstrahlung photon input will be negligible compared 
to cyclotron/synchrotron and the external soft photon fields. Thus, we assume u s = ii cx t + u sy 
with 

Uext = (A7) 

and 

1 U J A AT? °f 

n sy = - / duB u (e e ) + — J dvj v {@ e ). (A8) 

^SSA 

The photon energy I s is e ex t if ^ext > u sy , or maxjessA, eo} else, where eo = (hv c /m e c 2 ). 
The bremsstrahlung cooling rate is 

fdW\ = _( ^ar 2 e m e c 3 n e V@~ e for 6 e < 1 

V dt ) br ~ 1 96 a r\ m e c 3 n e 6 e (ln[26 e ] + 0.673) for G e ;> 1 



(Haug 1985) where a = 1/137 is the fine structure constant and r e is the classical electron radius. 



The heating rate due to Alfven wave turbulence is either specified as a fixed parameter, or 
can be evaluated using the acceleration rates according to Eq. (||) as 



fdW\ 
KoTJa 



oo 

2 



m e c z djjAfeh)- (A10) 



i 



Finally, the Coulomb heating / cooling rate is 
/ r1W\ l,m 

-7T = -— In An p ca T (feT p - kT e ) h(& e , 6 P ), (All) 
V ot J coui 1 m p 

where 

,,,„ Q , Vb; 2(e. + e p ) 2 + 2(e e + e p ) + i / i\ 

" (e " e " ) " ( e. + e,)^ IJTT) exp {-Wj (A12) 



( Dormer 1986[ ), In A is the Coulomb logarithm, n p is the number density of ions, and 



G p = kT p /(m p c 2 ), . 

The total heating/cooling rate consists the sum of all elementary processes, 

fdW\ __ /dW\ + /dW\ + /dW\ + /dW\ + /dW\ 
V dt / total V dt / Coul V dt J A V dt J sy V dt J 1C V dt / br ' 

and we are iteratively solving for the equilibrium solution(s) with (^-) = 0. 
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Fig. 1. — Equilibrium electron spectrum (a) and absolute value of the energy exchange coefficients 
(b) for our model calculation with parameters similar to the Cyg X-l case of Li et al. (1996b). 
The labels '+' and '-' in panel (b) indicate the sign of the respective contributions. The 
individual contributions are: synchrotron (thin solid), bremsstrahlung (thin long-dashed), Compton 
(thick long-dashed), Coulomb scattering (dot-dashed), M0ller + Bhabha scattering (short-dashed), 
hydromagnetic acceleration (dotted), and total (thick solid). Input parameters were: fcTj = 
3.5 MeV, r p = 0.7, R = 1.2- 10 8 cm, B = B cp = 1.1 *10 6 G, 5 2 = 0.059, kT BB = 0.1 keV, l s = 0.315. 
Resulting equilibrium parameters are I = 4.43, l st /l = 0.46, / pa i r = 1.7 %, kT e = 105 keV. The 
long-dashed curve in panel (a) is a thermal electron spectrum with kT e = 105 keV. 
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Fig. 2. — Equilibrium electron spectrum (a) and absolute value of the energy exchange coefficients 
(b) for our model calculation with parameters similar to the GRO J0422+32 case of Li et al. 
(1996b). The labels '+' and '-' in panel (b) indicate the sign of the respective contributions. 
The individual contributions are: synchrotron (thin solid), bremsstrahlung (thin long-dashed), 
Compton (thick long-dashed), Coulomb scattering (dot-dashed), M0ller + Bhabha scattering 
(short-dashed), hydromagnetic acceleration (dotted), and total (thick solid). Input parameters 
were: kT { = 7.2 MeV, r p = 0.7, i? = 1.2 ■ 10 8 cm, B = B cp = 1.65 * 10 6 G, 5 2 = 0.065, 
kT-BB = 0.1 keV, l s = 0.72. Resulting equilibrium parameters are I = 8.66, l s t/l = 0.43, / pa i r = 11 %, 
kT e = 93 keV. The long-dashed curve in panel (a) is a thermal electron spectrum with kT e = 93 keV. 
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Fig. 3. — Equilibrium electron spectra (top panel) and photon spectra (bottom panel) for fixed 
pulsar dipole moment ^30 = 1 and no plasma wave turbulence (5 2 = 0). 
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Fig. 5. — Equilibrium electron spectra (left panels) and photon spectra (right panels) for fixed 
pulsar dipole moment ^30 = 1. The legend in the upper left panel refers to the values of U used in 
the individual runs. 
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Fig. 6. — Same es Fig. ^, but with /U30 = 10 
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Fig. 7. — Equilibrium electron spectra (left panels) and photon spectra (right panels) for fixed 
accretion rate and luminosity I* = 0.25. The legend in the upper left panel refers to the values of 
/i30 used in the individual runs. 



